\documentclass[11pt]{article}
%\documentclass[11pt,draft]{article} 
%\pdfoutput=1 
%\usepackage{showkeys}
\usepackage[utf8x]{inputenc}
\usepackage[french]{babel}
%\usepackage[english]{babel}
\usepackage{graphicx}           
%\usepackage{url}
%\usepackage{eurosym}
%\usepackage{makeidx}
\usepackage[plain]{algorithm}
\usepackage{algorithmic} 
%\usepackage{color} 
%\usepackage[margin=3.3cm]{geometry}
%\usepackage{setspace}
%\usepackage{datetime}
%\usepackage{pifont}
%\usepackage{tikz}
\usepackage[numbers]{natbib}
%\usepackage{endfloat}

\usepackage{amsmath,amsfonts,amssymb,latexsym}
%stmaryrd}

%\doublespacing



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%---------------------------------------------------------------------
\input lib/f-macros
\input lib/macros
\input lib/f-algorithmic
\newenvironment{psmallmatrix}{\left(\begin{smallmatrix}}{\end{smallmatrix}\right)}
\newcommand{\Sin}{S_{\textrm{\tiny in}}}
\newcommand{\mumax}{\mu_{\textrm{\tiny max}}}
\newcommand{\Mfilter}{M_{\textrm{\tiny filtre}}}
\newcommand{\deltafilter}{\delta_{\textrm{\tiny filtre}}}
\newcommand{\deltasimu}{\delta_{\textrm{\tiny simu}}}
\newcommand{\Msimu}{M_{\textrm{\tiny simu}}}
\newcommand{\nmax}{n_{\textrm{\tiny max}}}
%---------------------------------------------------------------------


%---------------------------------------------------------------------
\title{Filtrage particulaire pour le chemostat}
\author{Fabien, Marc, Jérôme, Boumédiène etc.}
\date{{\small Créé: Juillet 29, 2011}\\{\small Dernière compilation: \today}}
%---------------------------------------------------------------------


%---------------------------------------------------------------------
\maketitle
%---------------------------------------------------------------------


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\begin{abstract}
%xxxxx
%\paragraph{Keywords and phrases:} 
%xxxx
%\end{abstract}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


Je présente la version la plus simple du filtre particulaire appliqué au plus simple chemostat: je suppose tous les paramètres connues. Il y a plusieurs pistes de travail:
\begin{enumerate}
\item faire ce filtrage en plus grande dimension;
\item estimation de paramètres (sur ce modèle ou un plus complexe);
\item mismatch: simuler avec un modèle, filtrer avec un autre, et être capable de dire quelque chose;
\item détection de changement: des caractéristiques de l'expérimentation on brusquement changé, comment le détecter, l'identifier etc. Cela sert aussi en exploitation industrielle.
\end{enumerate}
Pour la suite immédiate, je pense que Boumédiène devrait regarder \fenumi.

\smallskip

Le filtre a un intérêt majeur lorsque que l'on s'intéresse à ce qui se passe à l'instant courant et lorsque la fréquence des observations est relativement élevée.
Lorsque la fréquence est trop faible, il faudra passer la main à d'autres méthodes. Lorsque l'on s'intéresse à ce qui se passe le long de tout l'intervalle de temps, il serait intéressant de faire appel à des techniques de ``lissage'' ou à d'autres techniques.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Le problème}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

On considère la solution $Z_{t}=(X_{t},S_{t})$ du modèle stochastique du chemostat \cite{campillo2011chemostat} sous forme d'équation différentielle stochastique (EDS):
\begin{subequations}
\label{eq.state}
\begin{align}
\label{eq.state.X}
  \rmd X_{t} &= (\mu(S_{t})-D)\,X_{t}\,\rmd t
  +c_1\,\sqrt{X_{t}}\,\rmd W^1_{t}
\\
\label{eq.state.S}
  \rmd S_{t} &= -k\,\mu(S_{t})\,X_{t}\,\rmd t+D\,(\Sin-S_{t})\,\rmd t
  +c_2\,\sqrt{S_{t}}\,\rmd W^2_{t}
\end{align}
\end{subequations}
pour $t\in[0,T]$, où $X_{t}$ et $S_{t}$ désignent les concentrations (g/l) en biomasse et en substrat à l'instant $t$; $W_{t}^1$ et $W_{t}^2$ sont des mouvements browniens standard scalaires. Les paramètres sont: le taux de dilution  $D$ (1/h), la concentration du substrat en entrée $\Sin$ (mg/l), le coefficient st\oe chiométrique  
$k$ et les intensités  $c_{1}$ et $c_{2}$ des bruits. 

La densité de la loi de la condition initiale $Z_{0}=(X_{0},S_{0})$ est notée $p_{0}(z)=p_{0}(x,s)p_{Z_{0}}(z)=p_{X_{0},S_{0}}(x,s)$.


La fonction de croissance spécifique  
$\mu(s)$ est de type Monod ou Haldane:
\begin{align}
  \mu(s) &= \frac{\mumax\,s}{k_{s}+s} \quad\textrm{(Monod)}\,,
  &
  \mu(s) &= \frac{\mu_{0}\,s}{k_{s}+s+s^2/k_{i}} \quad\textrm{(Haldane)}\,.
\end{align}
Pour le modèle de Monod, les paramètres sont le taux de croissance maximum  $\mumax$ (1/h) et la constante de demi-saturation $k_{s}$ (mg/l). 

\bigskip

Nous disposons d'observations aux instants $t_{n}=n\,\Delta$:
\begin{align}
  Y_{n} &= S_{t_{n}} + \sigma\,S_{t_{n}} \,v_n
\end{align}
où $v_{n}\simiid N(0,1)$ et $T=\nmax\,\Delta$. Les bruits $W^i_{t}$ des équations d'état, les bruits d'observation $v_{n}$ et la condition initiale   $(X_{0},B_{0})$ sont supposés indépendants.

Nous utiliserons les notations 
$f_{1}(z)=f_{1}(x,s)=(\mu(s)-D)\,x$, $g_1(x)=c_1\,\sqrt{x}$,
 $f_{2}(x,s)=-k\,\mu(s)\,x+D\,(\Sin-s)$, et $g_2(s)=c_2\,\sqrt{s}$.

\medskip

À l'instant $t_{n}$, on dispose d'observations $y_{0},\dots,y_{n}$ (noté $y_{0:n}$) qui sont des réalisations de $Y_{0},\dots,Y_{n}$ (noté $Y_{0:n}$). Le problème de filtrage non-linéaire consiste à calculer le filtre non-linéaire optimal $\pi_{n}(z)=\pi_{n}(x,s)$ défini par:
\begin{align}
  \pi_{n}(z) \eqdef p_{Z_{t_{n}}|Y_{0:n}}(z|y_{0:n})\,.
\end{align}
Le filtre est donc la densité de $Z_{t_{n}}$ sachant les observations $Y_{0:n}=y_{0:n}$. Ce filtre permet par exemple de calculer l'estimation  $\hat Z_{t_{n}}$ de $Z_{t_{n}}$
\[
  \hat Z_{t_{n}}
  =
  \E(Z_{t_{n}}|Y_{0:n}=y_{0:n})
  =
  \int_{[0,\infty[^2} z\,\pi_{n}(z)\,\rmd z
\]
ou encore $(\hat X_{t_{n}},\hat S_{t_{n}})
= \int_{0}^\infty\int_{0}^\infty (x,s)\,\pi_{n}(x,s)\,\rmd x\,\rmd s$ ainsi que des estimations de la forme
$  \E(\phi(Z_{t_{n}})|Y_{0:n}=y_{0:n})
  =
  \int_{[0,\infty[^2} \phi(z)\,\pi_{n}(z)\,\rmd z
$, 
$  \P(Z_{t_{n}}\in A|Y_{0:n}=y_{0:n})
  =
  \int_{A} \phi(z)\,\pi_{n}(z)\,\rmd z
$, etc.

\medskip

On définit également le filtre prédit
\begin{align}
  \pi_{n^-}(z) 
  \eqdef
  p_{Z_{t_{n}}|Y_{0:n-1}}(z|y_{0:n-1})\,.
\end{align}
Le filtre prédit $\pi_{n^-}$ contient toute l'information disponible sur $Z_{t_{n}}$ à partir des observations $Y_{0:n-1}=y_{0:n-1}$, c'est-à-dire avant l'observation de $Y_n=y_{n}$ ; le filtre 
$\pi_{n}$ contient toute l'information disponible sur $Z_{t_{n}}$ avec toutes les observations $Y_{0:n}=y_{0:n}$.

\medskip


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Le filtre optimal}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Le filtre optimal va faire appel à deux outils: d'une part le \emph{noyau de transition} de l'équation d'état \eqref{eq.state} et d'autre part la \emph{fonction de vraisemblance} correspondant à la dernière observation; nous introduisons ces outils avant de décrire l'itération $\pi_{n-1} \to \pi_{n}$.

%---------------------------------
\subsection*{\it Noyau de transition}


On définit le noyau de transition de l'équation d'état \eqref{eq.state}:
\begin{align}
\label{eq.Q}
  Q_{\Delta}(z,z')
  \eqdef
  p_{Z_{t+\Delta}|Z_{t}}(z'|z)
  \,,
\end{align}
$z'\to Q_{\Delta}(z,z')$ est donc la densité de $Z_{t+\Delta}$ sachant que $Z_{t}=z$. Par exemple:
\[
   p_{Z_{t+\Delta}}(z')
   =
   \int_{[0,\infty[^2} Q_{\Delta}(z,z')\,p_{Z_{t}}(z)\,\rmd z
   =
   \int_{[0,\infty[^2} p_{Z_{t+\Delta}|Z_{t}}(z'|z)\,p_{Z_{t}}(z)\,\rmd z
\]
i.e. la noyau $Q_{\Delta}(z,z')$ ``transporte'' la densité  $p_{Z_{t}}$ pour la ``transformer'' en $p_{Z_{t+\delta}}$.

%---------------------------------
\subsection*{\it Fonction de vraisemblance}

On définit la fonction de vraisemblance locale $\Psi(y,z)$. Il s'agit de la  densité conditionnelle de l'observation $Y_{n}=y$ sachant que $Z_{t}=z=(x,s)$, c'est-à-dire:
\begin{align}
\label{eq.likelihood}
  \Psi(y,z)
  &\eqdef
  \exp\Big(-\frac{1}{2\,\sigma^2\,s^2}\,|y-s|^2\Big)\,.
\end{align}
Donc $\Psi(y,z)$ sera grand lorsque l'observation $y$ est \emph{vraisemblablement} issu de l'état $z$ et petit dans le cas contraire.

\emph{Le filtre prédit est donc le filtre à l'itération précédente transportée par le noyau de l'équation d'état.}

%---------------------------------
\subsection*{\it Itération $\pi_{n-1} \to \pi_{n}$}

Partant de $\pi_{n-1}$, le filtre $\pi_{n}$ est obtenu séquentiellement à l'aide des deux étapes classiques: 
\begin{enumerate}
\item Prédiction ---
On calcule d'abord le filtre prédit  $\pi_{n^-}$ en transportant le filtre par le noyau $Q_{\Delta}$, c'est-à-dire:
\begin{subequations}
\label{eq.nlf}
\begin{align}
\label{eq.nlf.prediction}
  \pi_{n^-}(z')
  &=
  \int_{[0,\infty[^2} Q_{\Delta}(z,z')\,\pi_{n-1}(z)\,\rmd z\,.
\end{align}
\item Correction ---
Afin de calculer le filtre $\pi_{n}$, on utilise la nouvelle observation $Y_{n}=y_{n}$ en posant:
\begin{align}
\label{eq.nlf.correction}
\pi_{n}(z')
  =
  \frac{\Psi(y_{n},z)\,\pi_{n^-}(z)}
        {\int \Psi(y_{n},z')\,\pi_{n^-}(z')\,\rmd z'}
\end{align}
\end{subequations}
que l'on note plus simplement:
\[
\textstyle
  \pi_{n}(z')
  \propto
  \Psi(y_{n},z)\,\pi_{n^-}(z)\,.
\]
\end{enumerate}
\emph{Dans l'étape de prédiction le filtre est transporté par le noyau donnant le filtre prédit, dans l'étape de correction ce filtre prédit est repondéré par la vraisemblance de la nouvelle observation.}

%---------------------------------
\subsection*{\it Initialisation}

À l'instant $t=0$, $\pi_{0^-}=\pi^0$ et
\begin{align}
\label{eq.nlf.correction.ini}
  \pi_{0}(z')
  &\propto
  \Psi(y_0,z)\,\pi_{0^-}(z)\,.
\end{align}

\bigskip

On dispose donc d'une représentation analytique exacte du filtre optimal, mais aucune des deux étapes \eqref{eq.nlf.prediction} et \eqref{eq.nlf.correction} ne peut être résolue explicitement (sauf dans le cas linéaire/gaussien et dans le cas de quelques exemples très particuliers). L'étape \eqref{eq.nlf.prediction} est linéaire mais l'étape  \eqref{eq.nlf.correction} est fondamentalement non-linéaire.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Le filtre particulaire}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

On introduit l'approximation  particulaire $\pi_{n}^N$ la plus simple du filtre non linéaire $\pi_{n}$, c'est-à-dire le ``bootstrap filter'' \cite{doucet2001a,campillo2006d}. L'idée est de proposer une approximation du filtre optimal sous la forme:
\[
 \pi_{n}(z) 
 \simeq 
 \pi_{n}^N(z) 
 \eqdef 
 \sum_{i=1}^N\delta_{\xi^i_{n}}(z)
\]
où $\xi^1_{n},\dots,\xi^N_{n}$ sont les particules et où $\delta_{\xi}(z)$ désigne la mesure de Dirac en $\xi$ (\footnote{Pour toute fonction $\phi$, $\int \phi(z)\,\delta_{\xi}(z)\,\rmd z=\phi(\xi)$.}) 
Le filtre particulaire est une méthode de Monte Carlo \emph{sequentielle} où le calcul de $\pi_{n}^N$ nécessite seulement $\pi_{n-1}^N$ et la nouvelle mesure  $Y_{n}=y_{n}$. 

Partant des particules $(\xi^i_{n-1})_{i=1\cdots N}$, on calcule les nouvelles particules  $(\xi^i_{n})_{i=1\cdots N}$ en deux étapes: prédiction (également appelée mutation) et
 correction (également appelée sélection).

%--------------------------------------------------------------------------------
\subsection*{\it Prédiction (mutation)}
%--------------------------------------------------------------------------------

L'opération de transport par le noyau $Q_{\Delta}$ est complexe, on la remplace ici par sa version \emph{empirique}: on se contente d'\emph{échantillonner} selon $Q_{\Delta}$. Cette opération, notée $\xi \sim Q_{\Delta}(z,\cdot)$, signifie que $\xi=Z_{\Delta}$ où $(Z_t)_{0\leq t\leq \Delta}$ est la solution de \eqref{eq.state} issue de $Z_0=z$ à l'instant $t=0$. 

On aura également besoin de ``l'échantillonneur de trajectoires'' $\QQ_{\Delta}(z,\zeta)$ où  $\zeta=(\zeta_{t})_{t\in[0,\Delta]}\sim \QQ_{\Delta}(z,\cdot)$ signifie que $(\zeta_{t})_{t\in[0,\Delta]}$ est la solution de   \eqref{eq.state} pour $t\in[0,\Delta]$ partant de $\zeta_{0}=z$ à l'instant $t=0$.

\medskip

Dans une première étape de prédiction, les particules $\xi^i_{n-1}$ ``singent'' indépendamment les unes des autres l'évolution du système  \eqref{eq.state}, i.e. on les fait évoluer selon le noyau de transition  $Q_{\Delta}$ défini par  \eqref{eq.Q}. Cela ne peut se faire qu'à l'aide de simulations numériques.  En pratique, on utilise un shéma d'Euler-Maruyama conduisant à un noyau approché  $\tilde Q_{\Delta}$, voir l'Algorithme \ref{algo.Q}.
Ainsi on échantillonne selon:
\[
   \tilde\xi^i_{n} \sim \tilde Q_{\Delta}(\xi^i_{n-1},\cdot)\,.
\]




%--------------------
\begin{algorithm}
\begin{center}
\begin{minipage}{\textwidth}
\hrulefill\\[-1em]
%\mbox{}
\begin{algorithmic}
\STATE $(x,s)$ donné
\STATE $\deltafilter \ot \Delta/\Mfilter$
\FOR {$m=1:\Mfilter$}
  \STATE $w_{1}\sim N(0,1)$, $w_{2}\sim N(0,1)$
  \STATE $x' \ot \max(0,x+f_{1}(x,s)\,\deltafilter+g_{1}(x)\,\sqrt{\deltafilter}\,w_{1})$
  \STATE $s' \ot \max(0,s+f_{2}(x,s)\,\deltafilter+g_{2}(s)\,\sqrt{\deltafilter}\,w_{2})$
  \STATE $x \ot x'$, $s \ot s'$
\ENDFOR
\end{algorithmic}
\hrulefill
\end{minipage}
\end{center}
\caption{\it Simulation de l'équation \eqref{eq.state} à l'aide d'un schéma d'Euler-Maruyama donnant une approximation $\tilde Q_{\Delta}$ du noyau  $Q_{\Delta}$ défini par \eqref{eq.Q}. Afin d'éviter les valeurs négatives de concentration, on utilise l'opérateur $\max(0,\cdot)$.}
\label{algo.Q}
\end{algorithm}
%--------------------



%----------------------------------------------------------------------------------
\subsection*{\it Correction (sélection)}
%----------------------------------------------------------------------------------

Dans une seconde étape de correction, à chacune des particules prédites on associe un poids proportionnel à la vraisemblance de la nouvelle observation, c'est-à-dire:
\[
   \omega^i_{n}
   \propto
   \Psi(y_{n},\xi^{i}_{n})
\]
avec $\sum_{i=1}^N\omega^i_{n}=1$. Les particules sont alors sélectionnées selon ses poids, en d'autres termes on échantillonne $N$ particules $\xi^1_{n},\dots,\xi^N_{n}$ selon la loi empirique $\sum_{i=1}^N\omega^i\,\delta_{\tilde\xi^i_{n}}$:
\begin{align*}
  \xi^1_{n},\dots,\xi^N_{n}
  \simiid
  \sum_{i=1}^N\omega^i_{n}\,\delta_{\tilde\xi^i_{n}}\,.
\end{align*}

\bigskip


\emph{En résumé, dans l'étape de prédiction les particules changent d'état, elles {\bfseries mutent},
en ``singeant'' l'équation d'état; dans l'étape de correction, on associe des poids de vraisemblance à chacune des particules, on effectue alors une {\bfseries sélection} favorisant les particules très vraisemblables au détriment des peu vraisemblable. Ce poids de vraisemblance est grand lorsque la nouvelle mesure $y_{n}$ aurait pu vraisemblament être issue de la particule associée.
} Cette version du filtre, appelée ``bootstrap'', est la plus simple et la plus naturelle, il existe de très nombreuses versions.

L'algorithme de filtrage particulaire est détaillé dans l'Algorithme  
\ref{algo.particle}.



%--------------------
\begin{algorithm}
\begin{center}
\begin{minipage}{\textwidth}
\hrulefill\\[-1em]
%\mbox{}
\begin{algorithmic}
\STATE{ }\COMMENT{initialisation}     
\STATE $\tilde\xi^{1},\dots,\tilde\xi^{N}\simiid \law(Z_{0})$
\STATE $\omega^i \ot \Psi(y_{0},\tilde\xi^{i})$
     pour $i=1:N$
     \COMMENT{poids de vraisemblance}
     \vphantom{$\int_{0}^N$}
\STATE $\omega^i \ot \omega^i/\sum_{i'=1}^N\omega^{i'}$
     pour $i=1:N$
     \COMMENT{normalisation}
     \vphantom{$\int_{0}^N$}
\STATE $\xi^{1},\dots,\xi^{N}\simiid \sum_{i=1}^N\omega^i\,\delta_{\xi^i}$
     \COMMENT{correction (sélection)}
     \vphantom{$\int_{0}^N$}
\STATE sauvegarder $(0,\xi^{1:N})$
     \vphantom{$\int_{0}^N$}
\STATE{ }     
\STATE{ }\COMMENT{itérations}     
\FOR {$n=1,2,\dots$}
  \STATE  $\tilde\xi^i \sim \tilde Q_{\Delta}(\xi^i,\rmd z')$
     pour $i=1:N$
     \COMMENT{prédiction (mutation)}
     \vphantom{$\int_{0}^N$}
  \STATE $\omega^i \ot \Psi(y_{n},\tilde\xi^{i})$
     pour $i=1:N$
     \COMMENT{poids de vraisemblance}
     \vphantom{$\int_{0}^N$}
  \STATE $\omega^i \ot \omega^i/\sum_{i'=1}^N\omega^{i'}$
     pour $i=1:N$
     \COMMENT{normalisation}
     \vphantom{$\int_{0}^N$}
  \STATE $\xi^{1},\dots,\xi^{N}\simiid \sum_{i=1}^N\omega^i\,\delta_{\xi^i}$
     \COMMENT{correction (selection)}
     \vphantom{$\int_{0}^N$}
  \STATE sauvegarder $(n\,\Delta,\xi^{1:N})$
     \vphantom{$\int_{0}^N$}
\ENDFOR
\end{algorithmic}
\hrulefill
\end{minipage}
\end{center}
\caption{\it Filtre particulaire (filtre ``bootstrap''): permet de calculer 
une approximation de $\pi_{n^-}$ de la forme $\sum_{i=1}^N\delta_{\tilde\xi^i}(\rmd z)$ et une approximation de $\pi_{n}(\rmd z)$ de la forme $\sum_{i=1}^N\delta_{\xi^i}(\rmd z)$.  Dans les étapes de prédiction, les particules prédites $\tilde\xi^i$ sont échantillonnés indépendamment les uns des autres.}
\label{algo.particle}
\end{algorithm}
%--------------------


\begin{remark}[correction progressive]
Dans l'Algorithme \ref{algo.particle}, lors\-que la somme des poids $\sum_{i'=1}^N\omega^{i'}$ est nulle cela signifie que le filtre a ``perdu la piste'' de l'état courant. 
Dans ce cas il existe plusieurs possibilités. En effet, dans la plupart des applications biotechnologiques la fréquence des observations est plutôt faible et il est donc possible de simuler des particules prédites jusqu'à obtenir un ensemble de particules ``suffisamment vraisemblables''. On peut également faire appel à des techniques spécifiques comme celle de la correction progressive, voir 
\citet[\S\,5.5]{oudjane2000b} et \citet{musso2001a}.
\end{remark}


\begin{remark}[Lissage]
Si  $T=\nmax\,\Delta$ est le temps final de l'expérimen\-tation ou de l'exploitation du bioréacteur, 
il peut être intéressant de calculer le lisseur: 
\[
  \bar\pi_{n}(z)
  \eqdef
  p_{Z_{n}|Y_{0:\nmax}}(z|y_{0:\nmax})\,,
\]
voir \cite{doucet2009tutorial,briers2010smoothing}
\end{remark}









%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Test numérique}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

On propose deux procédures \texttt{Matlab}:
\texttt{run\_simple()} dans \texttt{run\_simple.m},
et
\texttt{run()} dans \texttt{run.m}. Dans \texttt{run\_simple()}, 
le pas de temps des observations est égal au pas de temps de la simulation et à celui du filtre
 (i.e. $\Delta=\deltasimu=\deltafilter$ et $\Msimu=\Mfilter=1$). 
 Cela suppose que la fréquence des observations est grande comme dans le Test 1.
 Lorsque cette fréquence est basse, il est nécessaire de faire des itérations
 locales pour la simulation comme pour le filtre.
 
 \texttt{run\_simple()} ne génère pas de sorties. \texttt{run()} possède différents modes d'exécution:
\begin{description}
\item[\texttt{run(1)}] 
c'est le mode par défaut (i.e. correspondant à  \texttt{run()}); il exécute une simulation, le filtrage et un graphique temporel;
dans ce mode les données ne sont pas sauvées, c'est un mode d'exécution ``test''.
\item[\texttt{run(2)}] 
identique à   \texttt{run(1)} mais les données sont sauvegardées dans le fichier de données  
\verb+save_data.m+ contenant toutes les variables du test.
\item[\texttt{run(3)}] 
c'est le mode ``batch'' (sans graphique) qui peut être utilisé en mode batch sur calculateur, il génère le fichier de données \verb+save_data.m+.
\item[\texttt{run(4)}]
génère uniquement des graphiques (simulation, filtre), il fait appel au fichier  \verb+save_data.m+ dans le dossier courant afin de générer des graphiques:
\begin{verbatim}
   1. temporal plot t -> b(t) & s(t) 
   2. temporal plot t -> b(t) & s(t) with tubes 
   3. phase plot    (b(t),s(t)) for t in [0,T]
\end{verbatim}
(animations à venir). Avec \texttt{Matlab 7}, des fichiers \texttt{fig} et \texttt{pdf} sont générés   (\verb+graph_b_and_s+, \verb+graph_b_and_s_tubes+, \verb+graph_phase+), pour les versions plus anciennes des fichiers   \texttt{eps} sont générés.
\item[\texttt{run($i$)}] $i\neq 1,2,3,4$ est le mode ``aide''. 
\end{description}

%----------------------------------------------------------------------------------
\subsection{Simulation des observations}
%----------------------------------------------------------------------------------

Afin de filtrer, il est d'abord nécessaire de simuler des observations  $y_n$ à l'aide de l'Algorithme   \ref{algo.simu}.
Entre deux observations  $y_{n}$ on fait $\Msimu$ itérations locales.


%--------------------
\begin{algorithm}
\begin{center}
\begin{minipage}{\textwidth}
\hrulefill\\[-1em]
%\mbox{}
\begin{algorithmic}
\STATE $n \ot 0$, $t\ot 0$
\STATE $(x,s)\sim \law(Z_{0})$
\STATE $v\sim N(0,1)$
\STATE $y_{0}\ot s+\sigma\,s\,v$
\STATE $\deltasimu \ot \Delta/\Msimu$
\STATE sauver $(0,x,s)$ et $(0,y_{0})$
\FOR {$n=1:\nmax$}
\FOR {$m=1:\Msimu$}

  \STATE $t \ot t+\deltasimu$
  \STATE $w_{1}\sim N(0,1)$, $w_{2}\sim N(0,1)$
  \STATE $x' \ot \max(0,x+f_{1}(x,s)\,\deltasimu+g_{1}(x)\,\sqrt{\deltasimu}\,w_{1})$
  \STATE $s' \ot \max(0,s+f_{2}(x,s)\,\deltasimu+g_{2}(s)\,\sqrt{\deltasimu}\,w_{2})$
  \STATE $x \ot x'$, $s \ot s'$
  \STATE sauver $(t,x,s)$
\ENDFOR
\STATE $v\sim N(0,1)$
\STATE $y\ot s+\sigma\,s\,v$
\STATE sauver $(n\,\Delta,y)$
\ENDFOR

\end{algorithmic}
\hrulefill
\end{minipage}
\end{center}
\caption{\it Simulation du processus d'état $(X_{t},S_{t})$ et des observations   $Y_{n}=y_{n}$, $\Msimu$ est le nombre d'itérations locales pour la simulation de l'EDS entre deux instants d'observations; les sorties sont $(X_{t},S_{t})$ pour $t=n'\,\deltasimu$ et $n'=0,\dots, n\,\Msimu$, ainsi que  $(y_{n})_{n=0,\dots,\nmax}$.}
\label{algo.simu}
\end{algorithm}
%--------------------








%----------------------------------------------------------------------------------
\subsection{Test 1}
%----------------------------------------------------------------------------------


L'instant final est  $T=1000$ (h) et l'on utilise le même pas de discrétisation en temps pour la simulation de l'EDS, pour les observations et pour le filtre:
$\Delta=\deltasimu=\deltafilter=0.5$ (i.e. $\Msimu=\Mfilter=1$). Les paramètres sont: taux de dilution $D=0.01\,\textrm{h}^{-1}$; concentration du substrat en entrée  $\Sin=100\,\textrm{mg}/\textrm{l}$; coefficient st\oe chiométrique $k=10$; taux de croissance maximum  $\mumax=0.3\,\textrm{h}^{-1}$; coefficient de demi-saturation   $k_{s}=10\,\textrm{mg}/\textrm{l}$. L'instant final  $T=1000$ correspond à 10 fois le temps de rétention $\frac1D$.

Les intensités des bruits d'état sont $c_{1}=c_{2}=0.03$, l'intensité du bruit d'observation est $\sigma=0.2$. La loi initiale est  $\pi^0(x,s)=\NN(0.2,0.5^2)\otimes\NN(1,0.5^2)$.

Le nombre de particules est $1000$.




%-------------------
\begin{figure}
\centering
\includegraphics[width=11cm]{matlab/test1/graph_b_and_s.pdf}
\caption{Test 1, sorties obtenues avec  \texttt{run(4)}; le processus simulé    $t\to (X_{t},S_{t})$ (bleu), le processus d'observation  $t_{n}\to y_{n}$ (vert), les estimations   $t\to (\hat X_{t},\hat S_{t})$ (rouge) et la trajectoire déterministe  $t\to (x(t),s(t))$ (noir) (procédure graphique \texttt{graph\symbol{95}b\symbol{95}and\symbol{95}s}).}
\label{fig.example1}
\end{figure}
%-------------------

%-------------------
\begin{figure}
\centering
\includegraphics[width=11cm]{matlab/test1/graph_b_and_s_tubes.pdf}
\caption{Test 1, sorties obtenues à l'aide de \texttt{run(4)}; identique à la  Figure \ref{fig.example1} mais sans le processus d'observation et avec les ``tubes'' de confiance en gris autour des estimations  $\hat X_{t}$ et $\hat S_{t}$ correspondant aux valeurs minimum et maximum prisent par les composantes  $x$ et  $s$ des particules (procédure graphique \texttt{graph\symbol{95}b\symbol{95}and\symbol{95}s\symbol{95}tubes}).}
\label{fig.example2}
\end{figure}
%-------------------

%-------------------
\begin{figure}
\centering
\includegraphics[width=11cm]{matlab/test1/graph_phase.pdf}
\caption{Test 1, tracé dans l'espace des phases obtenu par  \texttt{run(4)} (procédure graphique \texttt{graph\symbol{95}phase}).}
\label{fig.example3}
\end{figure}
%-------------------


%----------------------------------------------------------------------------------
\subsection{Test 2}
%----------------------------------------------------------------------------------

Ce test est identique au Test 1 sauf que l'on n'utilise pas ici le même pas pour la simulation de l'EDS, pour les observations comme pour le filtre: $\Delta=10$, $\Delta=\deltasimu=\deltafilter=0.5$ (i.e. $\Msimu=\Mfilter=20$). 
Dans ce test, le filtre particulaire est nettement plus rapide que dans le Test 1, les temps CPU respectifs sont 17.3945 (s) pour le Test 1 et 0.84054 (s) pour le  Test 2. Cela est dû au fait que l'étape de sélection ne peut pas être ``vectorisée''
comme l'étape de prédiction, elle est donc gourmande en temps calcul.


%-------------------
\begin{figure}
\centering
\includegraphics[width=11cm]{matlab/test2/graph_b_and_s.pdf}
\caption{Test 2, sorties obtenues avec \texttt{run(4)}; processus simulé $t\to (X_{t},S_{t})$ (bleu), processus d'observation $t_{n}\to y_{n}$ (vert), les estimations  $t\to (\hat X_{t},\hat S_{t})$ (rouge) et la trajectoire déterministe  $t\to (x(t),s(t))$ (noir) (procédure graphique \texttt{graph\symbol{95}b\symbol{95}and\symbol{95}s}).}
\label{fig.test2.example1}
\end{figure}
%-------------------

%-------------------
\begin{figure}
\centering
\includegraphics[width=11cm]{matlab/test2/graph_b_and_s_tubes.pdf}
\caption{Test 2, sorties obtenues avec \texttt{run(4)}; identique à la  Figure \ref{fig.test2.example1} mais sans le processus d'observation et avec les ``tubes'' de confiance en gris autour des estimations  $\hat X_{t}$ et $\hat S_{t}$ correspondant aux valeurs minimum et maximum prisent par les composantes  $x$ et  $s$ des particules (procédure graphique \texttt{graph\symbol{95}b\symbol{95}and\symbol{95}s}).}
\label{fig.test2.example2}
\end{figure}
%-------------------

%-------------------
\begin{figure}
\centering
\includegraphics[width=11cm]{matlab/test2/graph_phase.pdf}
\caption{Test 2, tracé dans l'espace des phases obtenu par  \texttt{run(4)} (procédure graphique \texttt{graph\symbol{95}phase}).}
\label{fig.test2.example3}
\end{figure}
%-------------------


\clearpage
\appendix

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Notations de probabilités}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%---------------------------------------------------------------------
\subsection*{Variables aléatoires, (densités de) lois de variables aléatoires}
%---------------------------------------------------------------------

On introduit ici les notations probabilistes couramment utilisées en in\-génierie. Elles ne sont pas
mathématiquement parfaitement rigoureuse mais elles ont l'avantage d'être très simples à comprendre et à utiliser.

Une variable aléatoire (v.a.) est une grandeur $X$ qui dépend du résultat d'une expérience aléatoire. Elle peut être  à valeurs par exemple dans $\R$ et on parlera de v.a. réelle (v.a.r.); ou à valeurs dans $\R^n$ et on parlera de v.a. vectorielle.

La façon dont les valeurs que peut prendre une v.a. $X$ se distribuent dans l'espace $\R^n$ est caractérisée par une fonction  $x\mapsto p_{X}(x)$ appelée densité de la loi de $X$ (ou simplement densité de $X$ ou encore, par abus, loi de $X$). Il s'agit d'une fonction positive telle que $\int_{\R^n} p_{X}(x)\,\rmd x=1$ (i.e. qui somme à 1) et telle que:
\[
  \P(X\in A)
  =
  \int_{A}p_{X}(x)\,\rmd x 
\]
pour tout $A\subset\R^n$; cette loi permet également de calculer la valeur moyenne d'une fonction $\phi(X)$ de $X$ à l'aide de l'espérance $\E\phi(X)$ définie par:
\[
  \E(\phi(X))
  =
  \int_{\R^n}\,\phi(x)\,p_{X}(x)\,\rmd x 
\]
pour tout fonction $\phi:\R^n\mapsto\R$ (\footnote{Mathématiquement, il existe des lois qui n'admettent pas de densité, de plus $\P(X\in A)$ n'est défini que pour des $A$ qui sont des produits d'intervalles ou des réunions dénombrables de telles ensembles (encore appelés ensembles boréliens); de même $\E\phi(X)$ n'est définit que pour des fonctions $\phi$ mesurables et telles que $\int_{\R^n}\,|\phi(x)|\,p_{X}(x)\,\rmd x<\infty$, donc par exemple pour des fonctions mesurables et bornées.}). Noter que si on introduit la fonction indicatrice $\indic_{A}$ d'un ensemble $A$ définie par $\indic_{A}(x)=1$ si $x\in A$ et 0 sinon, alors $\E \indic_{A}(X)=\P(X\in A)$.

\medskip

Si $p(x,y)$ est la loi d'un couple $(X,Y)$ de v.a. à valeurs dans $\R^n\times\R^d$, alors la densité marginale de $X$ est définie par $p_{X}(x)=\int_{\R^d} p(x,y)\,\rmd y$. Il s'agit donc de la densité de $X$, cette fonction contient toute la connaissance que l'on a sur la v.a. de $X$. Supposons maintenant que l'on observe $Y=y$, i.e. que l'on observe une réalisation particulière de $Y$. Comment cette information peut-elle être utilisée pour améliorer notre connaissance de $X$~? La réponse est donnée par la (densité de la) loi conditionnelle de $X$ sachant $Y=y$, notée $p_{X|Y}(x|y)$ ou encore $p_{X|Y=y}(x)$, définie par:
\[
   p_{X|Y}(x|y)
   \eqdef
   \frac
   {p(x,y)}
   {p_{X}(x)}
=
   \frac
   {p(x,y)}
   {\int_{\R^d}p(x',y)\,\rmd x'}\,.
\] 
À tout $y$ fixé, il s'agit d'une densité puisque $x\to p_{X|Y}(x|y)$ est une fonction positive qui somme à 1.
La probabilité conditionnelle de $X\in A$ sachant que $Y=y$ est donnée par
\[
  \P(X\in A|Y=y)
  =
  \int_{A}p_{X|Y}(x|y)\,\rmd x 
\]
pour tout $A\subset\R^n$
et l'espérance conditionnelle de $\phi(X)$ sachant $Y=y$ est définie par
\[
  \E(\phi(X)|Y=y)
  =
  \int_{\R^n}\,\phi(x)\,p_{X|Y}(x|y)\,\rmd x 
  \,.
\]
On dit que les v.a. $X$ et $Y$ sont indépendantes, noté $X\indep Y$, lorsque $p(x,y)=p_{X}(x)\,p_{Y}(x)$ et alors: $\P(X\in A,Y\in B)=\P(X\in A)\,\P(Y\in B)$ et $\E(\phi(X)\,\psi(X))=\E(\phi(X))\,\E(\psi(X))$. Lorsque $X\indep Y$ alors $p_{X|Y}(x|y)=p_{X}(x)$, c'est-à-dire : l'observation de $Y=y$ n'apporte aucune information supplémentaire sur la connaissance que l'on a de $X$, d'où la terminologie ``indépendance''.

%---------------------------------------------------------------------
\subsection*{Variables discrètes et leur simulation}
%---------------------------------------------------------------------

[à faire]

%---------------------------------------------------------------------
\subsection*{Processus, loi de processus, processus de Markov}
%---------------------------------------------------------------------

Un processus aléatoire est une collection $(X_{t})_{t\in\TT}$ de v.a. indicées
par un ensemble qui peut être discret, $\TT=\N$, et on parlera de processus en temps discret, ou continu,
$\TT=[0,T]$ ou $[0,\infty[$, et on parlera de processus en temps continu.

\medskip

La loi d'un processus en temps discret $(X_n)_{n\in\N}$ est donné par la données
des lois $p_{X_{0},\dots,X_{n}}(x_{0},\dots,x_{n})$ des vecteurs $(X_{0},\dots,X_{n})$ pour tout $n$. La factorisation suivante est toujours valide:
\begin{align*}
&
p_{X_{0},\dots,X_{n}}(x_{0},\dots,x_{n})
=
\\
&\qquad
=
p_{X_{n}|X_{0},\dots,X_{n-1}}(x_{n}|x_{0},\dots,x_{n-1})
\\
&\qquad\qquad
\times p_{X_{n-1}|X_{0},\dots,X_{n-2}}(x_{n-1}|x_{0},\dots,x_{n-2})\times
\cdots
\\
&\qquad\qquad\qquad\qquad\qquad
\cdots\times
p_{X_{1}|X_{0}}(x_{1}|x_{0})\times
p_{X_{0}}(x_{0})
\end{align*}
Ce processus est appelé \emph{chaîne de Markov} lorsque:
\begin{align*}
p_{X_{n}|X_{0},\dots,X_{n-1}}(x_{n}|x_{0},\dots,x_{n-1})
=
p_{X_{n}|X_{n-1}}(x_{n}|x_{n-1})
\end{align*}
c'est-à-dire l'état futur du processus ne dépend du passé que par l'intermédiaire 
du présent. Dans ce cas, la factorisation précédente se réduit~à:
\begin{align*}
p_{X_{0},\dots,X_{n}}(x_{0},\dots,x_{n})
&
=
p_{X_{n}|X_{n-1}}(x_{n}|x_{n-1})
\\
&\qquad\qquad
\times p_{X_{n-1}|X_{n-2}}(x_{n-1}|x_{n-2})\times
\cdots
\\
&\qquad\qquad\qquad
\cdots\times
p_{X_{1}|X_{0}}(x_{1}|x_{0})\times
p_{X_{0}}(x_{0})
\end{align*}
Donc la loi du processus de Markov ne dépend que de $p_{X_{n}|X_{n-1}}(x_{n}|x_{n-1})$ appelée ``noyau de transition'' et de la loi initiale $p_{X_{0}}(x_{0})$.

\medskip

La loi  d'un processus en temps continu $(X_{t})_{t\geq 0}$ est entièrement caractérisée par la loi des vecteurs $(X_{t_{0}},X_{t_{1}},\dots,X_{t_{n}})$ pour tout $n$ et tout $0=t_{0}<t_{1}<\cdots < t_{n}$. À nouveau on a:
\begin{align*}
&
p_{X_{t_0},\dots,X_{t_n}}(x_{0},\dots,x_{n})
=
\\
&\qquad
=
p_{X_{t_n}|X_{t_0},\dots,X_{t_{n-1}}}(x_{n}|x_{0},\dots,x_{n-1})
\\
&\qquad\qquad
\times p_{X_{t_{n-1}}|X_{t_0},\dots,X_{t_{n-2}}}(x_{n-1}|x_{0},\dots,x_{n-2})\times
\cdots
\\
&\qquad\qquad\qquad\qquad\qquad
\cdots\times
p_{X_{t_1}|X_{t_0}}(x_{1}|x_{0})\times
p_{X_{t_0}}(x_{0})\,.
\end{align*}
Ce processus est dit de Markov lorsque:
\begin{align*}
p_{X_{t_n}|X_{t_0},\dots,X_{t_{n-1}}}(x_{n}|x_{0},\dots,x_{n-1})
=
p_{X_{t_n}|X_{t_{n-1}}}(x_{n}|x_{n-1})
\end{align*}
c'est-à-dire l'état futur du processus ne dépend du passé que par l'intermédiaire 
du présent. Dans ce cas, la factorisation précédente se réduit~à:
\begin{align*}
&
p_{X_{t_0},\dots,X_{t_n}}(x_{0},\dots,x_{n})
=
p_{X_{t_n}|X_{t_{n-1}}}(x_{n}|x_{n-1})
\\
&
\qquad\qquad\qquad
\times p_{X_{t_{n-1}}|X_{t_{n-2}}}(x_{n-1}|x_{n-2})\times
\cdots
\\
&\qquad\qquad\qquad\qquad\qquad\qquad\qquad
\cdots\times
p_{X_{t_1}|X_{t_0}}(x_{1}|x_{0})\times
p_{X_{t_0}}(x_{0})\,.
\end{align*}
ainsi la loi du processus est entièrement caractérisée par le noyau de transition $p_{X_t|X_s}(x'|x)$ pour tout $s<t$, et par la loi initiale $p_{X_{t_0}}$.

Un exemple de processus de Markov en temps continu est donné par le mouvement brownien standard réel $(W_{t})_{t\geq 0}$, encore appelé processus de Wiener. Sa loi est définie par $p_{X_{t_0}}=\delta_{0}$ (i.e. $W_{0}=0$) et par:
\[
p_{W_t|W_s}(x'|x) 
=
g_{\NN(x,t-s)}(x)
\]
(la densité gaussienne de moyenne $x$ et de variance $t-s$), plus précisément le mouvement est un processus nulle en 0; à accroissement indépendants, i.e. pour tout $0=t_{0}<t_{1}<\cdots < t_{n}$ on a $(W_{t_{n}}-W_{t_{n-1}}) \indep (W_{t_{n-1}}-W_{t_{n-2}}) \indep \cdots \indep (W_{t_{1}}-W_{t_{0}})$; à accroissements gaussiens tels que $W_{t}-W_{s}\sim \NN(0,t-s)$ pour tout $s<t$.


%---------------------------------------------------------------------
\subsection*{Processus de diffusion, équation différentielle stochastique}
%---------------------------------------------------------------------

Une équation différentielle stochastique
\begin{align}
\label{eq.eds}
  \rmd X_{t}
  =
  f(X_{t})\,\rmd t
  +
  g(X_{t})\,\rmd W_t
\end{align}
est une notation pour :
\[
  X_{t}
  =
  X_0
  +
  \int_{0}^t f(X_{s})\,\rmd s
  +
  \int_{0}^t g(X_{s})\,\rmd W_s
\]
la première de ces intégrales est l'intégrale de Lebesgue classique, la seconde est l'intégrale stochastique de Itô, c'est la limite suivante:
\[
    \int_{0}^t g(X_{s})\,\rmd W_s
   \eqdef
   \lim_{n\to\infty}
   \sum_{i=0}^{n-1} g(X_{t_{i}})\,(W_{t_{i+1}}-W_{t_{i}})
\]
(au sens de la limite quadratique) avec $t_{i}=t\,\frac{i}{n}$. Cela signifie que l'EDS \eqref{eq.eds} est la limite du schéma d'Euler suivant (ou que l'EDS est une idéalisation du schéma suivant):
\[
  X_{t+\delta}
  =
  X_{t}
  +
  f(X_{t})\,\delta
  +
  g(X_{t})\,(W_{t+\delta}-W_{t})
\]
lorsque $\delta \to 0$. Noter que $W_{t+\delta}-W_{t}\sim \NN(0,\delta)$ et est indépendant des autres accroissements de $W_{t}$.

\medskip

L'EDS est bien posée et admet une unique solution lorsque les fonctions $f$ et $g$ sont à croissance au plus linéaire et sont localement lipschtiziennes, voir \cite{allen2007a} pour plus de rigueur. La solution de cette EDS est un processus de diffusion de coefficient de dérive $f$ et de coefficient de diffusion $g$.

\medskip

Noter que:
\begin{align*}
 \E \int_{0}^t g(X_{s})\,\rmd W_{s}&=0\,,
 \\
 \E\Big[ \Big(\int_{0}^t g(X_{s})\,\rmd W_{s}\Big)^2\Big]&=\E\Big[ \int_{0}^t g^2(X_{s})\,\rmd s\Big]\,.
\end{align*}
En particulier d'après \eqref{eq.eds}, on a:
\[
  \E X_{t}
  =
  \E X_0
  +
  \int_{0}^t \E f(X_{s})\,\rmd s
\]
et \textbf{si} $\E f(X_{s})=f(\E (X_{s}))$ alors en posant $x(t)=\E(X_{t})$, on a $\dot x(t)=f(x(t))$ mais cela n'est essentiellement vrai que dans le cas linéaire où $f(x)=\alpha\,x+\beta$, le plus souvent la solution $x(t)$ de l'équation différentielle n'est pas la moyenne du processus $X_{t}$.


%%---------------------------------------------------------------------
%\subsection*{Formule d'Itô}
%%---------------------------------------------------------------------
%
%\begin{align*}
%  \rmd \phi(X_{t})
%  &=
%  \phi'(X_{t})\,\rmd X_{t}
%  +
%  \demi\,\phi''(X_{t})\,g^2(X_{t})\,\rmd t
%\\
%  &=
%  \big[\phi'(X_{t})\,f(X_{t})+\demi\,\phi''(X_{t})\,g^2(X_{t})\big]\,\rmd t
%  +
%  \phi'(X_{t})\,d(X_{t})\,\rmd W_{t}
%\end{align*}
%par rapport au calcul différentiel classique, ce calcul fait apparaître une terme supplémentaire de Itô
%$\demi\,\phi''(X_{t})\,g^2(X_{t})\,\rmd t$.

%---------------------------------------------------------------------
\subsection*{Simulation d'une EDS}
%---------------------------------------------------------------------

Il suffit d'utiliser la construction même de l'intégrale stochastique. Sur une grille d'instant $t_{n}=n\,\Delta$ avec $\Delta=T/\nmax$, et l'approximation
en temps discret $(\tilde X_{t_{n}})_{n=0:\nmax}$ du processus $(X_{t})_{0\leq t\leq T}$ est obtenu en posant:
\[
  \tilde X_{t_{n+1}}
  =
  \tilde X_{t_{n}}
  +
  f(\tilde X_{t_{n}})\,\delta
  +
  g(\tilde X_{t_{n}})\,\sqrt{\delta}\,w_{n}
\]
où $w_{n}$ est une suite i.i.d. $N(0,1)$ en dimension 1 ou $N(0,I)$ en dimension supérieur.

Ce schéma est appelé Euler-Maruyama car c'est un schéma d'Euler et la preuve de sa convergence est due à Gisiro Maruyama.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\nocite{venkateswarlu2004a}
\addcontentsline{toc}{section}{Reference}
\bibliographystyle{plainnat}
\bibliography{lib/fab}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



